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The understanding of spin dynamics and relaxation mechanisms in clean graphene and the upper time and 
length scales on which spin devices can operate are prerequisites to realizing graphene spintronic technolo¬ 
gies. Here we theoretically reveal the nature of fundamental spin relaxation mechanisms in clean graphene 
on different substrates with spin-orbit Rashba fields as low as a few tens of peV. Spin lifetimes ranging from 
50 picoseconds up to several nanoseconds are found to be dictated by substrate-induced electron-hole charac¬ 
teristics. A crossover in the spin relaxation mechanism from a Dyakonov-Perel type for Si02 substrates to a 
broadening-induced dephasing for hBN substrates is described. The energy dependence of spin lifetimes, their 
ratio for spins pointing out-of-plane and in-plane, and the scaling with disorder provide a global picture about 
spin dynamics and relaxation in ultraclean graphene in presence of electron-hole puddles. 

PACS numbers: 72.80.Vp, 73.63.-b, 73.22.Pr, 72.15.Lh, 61.48.Gh 


The tantalizing prospect of graphene spintronics was initi¬ 
ated by Tombros and coworkers IT], who first reported long 
spin diffusion length in large area graphene. The small spin- 
orbit coupling (SOC) in carbon, plus the absence of a hyper- 
fine interaction, suggested unprecedented spin lifetimes (Ts) 
at room temperature (from /is to ms) |l2l [3 ^ |5] [h] E). 

However, despite significant progress in improving 
graphene quality, resolving contact issues, and reducing sub¬ 
strate effects III] 11 nni [m El El m [la , the measured 
Ts are orders of magnitude shorter, even for high-mobility 
samples. Extrinsic sources of SOC, including adatoms 
ifThl fm [TSl [T^ or lattice deformations 1^121], have been 
proposed to explain this discrepancy. Moreover, the na¬ 
ture of the dominant spin relaxation mechanism in graphene 
is elusive and debated. The conventional Dyakonov-Perel 
(DP) ll^ and Elliot-Yafet (EY) ll^ mechanisms, usually 
describing semiconductors and disordered metals, remain 
inconclusive in graphene because neither effect can con¬ 
vincingly reproduce the observed scaling between Ts and 
the momentum relaxation time Tp 10 El- Although gen¬ 
eralizations of both mechanisms have been proposed, they 
do not allow an unambiguous interpretation of experiments 

It should be noted that the achieved room-temperature 
spin lifetime in graphene is already long enough for the ex¬ 
ploration of spin-dependent phenomena such as the spin Hall 
effect Il26ll27]l. or to harness proximity effects as induced for 
instance by magnetic oxides ll28l or semiconducting tung¬ 
sten disulphide ll29l . However, a comprehensive picture of 
spin dynamics of massless Dirac fermions in presence of 
weak spin-orbit coupling fields is of paramount importance 
for further exploitation and manipulation of spin, pseudospin 
and valley degrees of freedom lEl [31] |E1 ■ 

In this study, we show numerically that a weak uniform 
Rashba SOC (tens of /reV), induced by an electric field or 


the substrate, yields spin lifetimes from 50 ps up to sev¬ 
eral nanoseconds. The dominant spin relaxation mechanism 
is shown to be dictated by long range potential fluctuations 
(electron-hole puddles) 1^ . Eor graphene on a Si02 sub¬ 
strate, such disorder is strong enough to interrupt the spin 
precession driven by the uniform Rashba field, resulting in 
motional narrowing and the DP mechanism. We also find the 
ratio /t| ~ 1 /2, demonstrating the anisotropy of the in¬ 
plane Rashba SOC field. Eor the case of a hexagonal boron 
nitride (hBN) substrate, the role of electron-hole puddles is 
reduced to an effective energy broadening and the spin life¬ 
time is limited by pure dephasing ll3^l35]l . These situations, 
however, share a common fingerprint - a M-shape energy 
dependence of Tg that is minimal at the Dirac point. Taken 
together, our results provide deeper insight into the funda¬ 
mentals of spin lifetimes in graphene dominated by electron- 
hole puddles. 


Results 

Disorder and Spin dynamics. 

Electron-hole puddles are real-space fluctuations of the 
chemical potential, induced by the underlying substrate, 
which locally shift the Dirac point li33 [36l [32l. Since 
measured transport properties usually result from an aver¬ 
age around the charge neutrality point, it is generally dif¬ 
ficult to access the physics at the Dirac point. As shown 
by Adam and coworkers 1331 . electron-hole puddles can be 
modeled as a random distribution of long range scatterers, 
= E^ie7 exp[-(r - ^j^ere ^ = 10 

and 30 nm denote the effective puddle ranges for Si02 and 
hBN substrates, respectively l3^ IMll . and Cj is randomly 
chosen within [—A, A]. Based on experimental data, typical 
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impurity densities are rii = 10^^ cm“^ (0.04%) for Si02 
and rii = 10^^ cm“^ (0.004%) for hBN substrates lf3^[^ . 
From such information, we can tune A to obtain suitable 
disorder profiles for the onsite energy of the 7r-orbital. Fig. 

(main frame) shows the onsite energy distribution corre¬ 
sponding to hBN and Si02 substrates, where the Gaussian 
profiles give standard deviations of cr = 5.5 and 56 meV, re¬ 
spectively. This allows us to extract A = 50 meV for Si02 
and A = 5 meV for hBN. The inset of Fig. shows the en¬ 
ergy landscape for a sample with 0.04% Gaussian impurities 
(Si 02 case). 
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FIG. 2: (a) Transport times for graphene on Si02 and hBN sub¬ 
strates (solid black and red curves, respectively). The dashed line 
shows the spin precession time, (b) Time-dependent spin polariza¬ 
tion for out-of-plane (solid red line) and in-plane (solid black line) 
spin injection for the Si 02 substrate, plus the fits to the exponential 
damping (dashed lines). The blue curves show the same informa¬ 
tion for the hBN substrate with out-of-plane injection. 


density n by n* = |n| -F 2^2 lf33l l42l l43l . The computed 
Tp are shown in Fig. [^a) for both substrates. For Si02, Tp 
is on the order of a few ps, while for hBN Tp is more than 
two orders of magnitude larger. The spin precession time, 
Tfi = irh/Xfi, is shown for comparison. 


FIG. 1: Onsite energy distribution of the carbon atoms in the 
graphene sample, which mimics the chemical potential induced by 
hBN (green) and Si02 (black) substrates together with their Gaus¬ 
sian fitting lines. Inset: Real space vizualization of the energy 
landscape for a graphene sample with 0.04% Gaussian impurities 
(Si 02 case). 


To fully characterize the role of electron-hole puddles, we 
evaluate the transport time Tp using a real-space order-N ap¬ 
proach, which computes the diffusion coefficient D{E,t). 
We extract Tp from the saturation of D{E,t) since Tp = 
Dmax{E)/ 2 vp{E) II 40 II . For numerical convenience, the 
calculations are made using A = O.I70 (with 70 the nearest 
neighbor hopping parameter), thus in absence of intervalley 
scattering ED, while the final values for hBN and Si02 sub¬ 
strates are extrapolated from the numerical results using the 
scaling law 



^-Kn-e 


( 1 ) 


where Ii{x) is the modified Bessel function of the first 
kind, Kq = 40.5ni(A/f)^(^/v/3a)^ is a dimensionless pa¬ 
rameter dictating the strength of the Gaussian potential, and 
the carrier density n* is modified from the pristine graphene 


Spin dynamics and lifetimes in the presence of electron-hole 
puddles. 

We now analyze the spin dynamics for puddles corre¬ 
sponding to the Si02 and hBN substrates. The blue curve in 
Fig.2(b) shows the time-dependent spin polarization for the 
hBN substrate (rii = 0.004%) at the Dirac point for an initial 
out-of-plane polarization, (see Methods). The po¬ 

larization exhibits oscillations with period Tq = 'kTi/Xr ~ 
55 ps, corresponding to the spin precession induced by the 
Rashba field. Simultaneously, the polarization decays in 
time, and by fitting = cos (27rf/Tn) both 

Tq and the spin relaxation time Ts can be evaluated. 

Fig.2(b) also shows (f) for the Si02 substrate (m = 

0.04%) with initial spin polarization in-plane (a =||) and 
out-of-plane (a =-L). In contrast to the hBN case, for 
which Pp^^{t) exhibits significant precession, the disorder 
strength of electron-hole puddles for Si 02 is sufficient to in¬ 
terrupt spin precession. As a result, the polarization for Si02 
is better fit with {t) = . The absence of preces¬ 

sion for Pp^^{t) compared to Pp^^{t) is consistent with 
the ratio between transport time and precession frequency, 
since Tp^'^^ /Tq 1 whereas > 1. 

To scrutinize the origin of the dominant relaxation mecha¬ 
nism, we first examine the spin lifetimes Tg for the Si 02 case 
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FIG. 3: Spin lifetimes for out-of-plane (a) and in-plane (b) spin 
injection for Si02 substrate at impurity densities of 0.04% (black 
solid curves), 0.08% (red dashed curves), and 0.16% (blue dotted 
curves ). (c) Spin lifetime with out-of-plane spin injection for the 
hBN substrate at impurity densities of 0.004% (black curve) and 
0.016% (red curve). 


when rotating spin polarization (out-of-plane vs. in-plane), 
and varying the density of impurities (0.04%, 0.08%, and 
0.16%). Fig.3 shows the extracted in the out-of-plane (a) 
and in-plane (b) cases. The energy dependence of Tg ex¬ 
hibits an M-shape increasing from a minimum at the Dirac 
point, with a saturation and downturn of for E > 200 
meV. The values of Tg range from 50 to 400 ps depending on 
the initial polarization and impurity density. We observe an 
increase of Tg with rii, which shows that a larger scattering 
strength reduces spin precession and dephasing, resulting in 
a longer spin lifetime, as described by the so-called motional 
narrowing effect ll44l . Additionally, the ratio t^/t| (not 
shown) changes from 0.3 to 0.45 when is varied from 
0.04% to 0.16%. Such behavior is expected when enhanced 
scattering drives more randomization of the direction of the 
Rashba SOC field, which ultimately yields /t]' = 0.5 in 
the strong disorder limit Gia. These results are fully con¬ 
sistent with the DP spin relaxation mechanism ll^ 12111441 . 

Fig.3(c) shows for the hBN substrate {rii = 0.004% 
and 0.016%) where a similar M-shape is observed. While 
t^{BN) is similar to r^(Si02) near the Dirac point, it is 
much larger at higher energies, reaching nearly 1 ns (for 
Xu = 37.4 peV). A striking difference is that the scaling 
of Tg with Ui is opposite to that of the Si 02 case, with an in¬ 
crease in puddle density resulting in a decrease in Tg, which 
indicates a different physical origin. For hBN, this behav¬ 
ior is reminiscent of the EY mechanism, but we will argue 
below that its origin is a different one. 


Crossover in spin relaxation behavior for hBN and Si02 
substrates. 

Fig.4 provides a global view of our results, where we plot 
Tg VS. 1/rp for the Si02 and hBN substrates (black and 
red symbols respectively) at the Dirac point and E = 
—200 meV (closed and open symbols respectively). For 
low defect densities (hBN substrate), Tg decreases strongly 
with decreasing Tp. However, with increasing defect den¬ 
sity (Si 02 substrate) this trend reverses and Tg scales al¬ 
most linearly with l/rp, according to the DP relationship 
Tg = V ■ (Tn/27r)^/Tp. ki E = —200 meV, v = 1, fitting 
the usual DP theory. At the Dirac point, the scaling is some¬ 
what weaker, with = 1/4. These results are reminiscent of 
those summarized in Fig. 5(a) of Drogeler et al. m, where 
spin lifetimes of graphene devices on Si 02 scaled inversely 
with the mobility, while devices on hBN appear to show the 
opposite trend. 

While the Si02 results of Fig.4 show DP behavior, the 
nature of the spin relaxation for weak electron-hole pud¬ 
dles is less clear. The fact that Tg and Tp decrease together 
suggests the EY mechanism, but we find Tg < Tp near the 
Dirac point and Tg Tp at higher energies. This contrasts 
with the usual picture of EY relaxation, where charge car¬ 
riers flip their spin when scattering off impurities, giving 
Ts = Tp/a, where a <C 1 is the spin flip probability 0. 
Instead, this situation matches that described in Ref. Il44l : 
when Tp > Tq, the spin precesses freely until phase infor¬ 
mation is lost during a collision, in analogy to the collisional 
broadening of optical spectroscopy. More collisions result 
in a greater loss of phase, reducing Tg with decreasing Tp. 
We verify this by removing the real-space disorder (setting 
A = 0) and modeling the electron-hole puddles with an ef¬ 
fective Lorentzian energy broadening rj*. The results are 
shown in Eig.4 (main frame, blue dashed line), where we 
plot Tg vs. rj* at E = —200 meV (top axis). Eor small 
rf, the scaling matches well with the real-space simulations 
of hBN, indicating that the puddles can be represented as 
a uniform energy broadening (See supplementary material). 
Larger values of rj* lead to stronger mixing of different spin 
dynamics and Tg saturates at very large rj*. There, the scal¬ 
ing of Tg vs. rj* clearly fails to replicate the DP behavior seen 
in the real-space simulations, since the effective broadening 
model does not induce the momentum scattering necessary 
for motional narrowing Il44l . 

We finally explain the downturn of Tg at the high energy 
wings of the M-shaped Tg behavior in the hBN case. We 
compare the spin dynamics in the TB model (Eq. in 
Methods) and the low-energy model in the absence of pud¬ 
dles (A = 0). In this regime Tp Tq, and spin dephas¬ 
ing and relaxation are driven by a combination of energy 
broadening and a nonuniform spin precession frequency. Eor 
the TB model, spin dynamics are calculated with the real- 
space approach and with a standard /c-space approach and 
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FIG. 4: Low-energy spin lifetimes versus 1 /rp (for initial out-of- 
plane spin polarization). Squares (circles) are for graphene on hBN 
(Si02) substrate. Closed (open) symbols are for spin relaxation at 
the Dirac point (at E — —200 meV). The blue dashed line shows 
the spin lifetime assuming only energy broadening (top axis). Inset: 
spin lifetime in absence of puddles computed using the TB model 
in real space (red circles) or fc-space (blue solid line), and the low- 
energy model in fc-space (green dashed line), with rj = 13.5 meV. 
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give identical Ts (inset of Fig.4, red circles and blue solid 
line), indicating the equivalence of the real- and A:-space 
approaches in the clean limit when accounting for the full 
TB Hamiltonian. We observe that while for all models, the 
spin lifetime shows a minimum at the Dirac point (in agree¬ 
ment with experimental data, and explained by a strong spin- 
pseudospin coupling in mx spin transport simulations 
with the widely used low-energy Hamiltonian (see 

Methods for and green dashed line in Fig. 4 inset for re¬ 
sults) clearly cannot capture the saturation and downturn of 
Ts{E), i.e. its full M-shape. To qualitatively reproduce the 
M-shape of Tg{E), the first-order term of the Rashba Hamil¬ 
tonian, ^^[kx{(yxSy + (JySx) + ky{axSx — aySy)], needs to 
be included in This term introduces stronger dephas¬ 

ing at higher energy, driven by the anisotropy of the Rashba 
spin-orbit interaction ifTSll . 

In addition to their different energy dependence, the TB 
and low-energy models also yield very different spin life¬ 
times. A value of = 10 ns is obtained at the Dirac point 
for the low-energy model, which is two orders of magnitude 
larger than Tg from the TB Hamiltonian, indicating a strong 
spin dephasing induced by the high-order fc-terms. Inter¬ 
estingly, by studying the changes of Tg{E) with respect to 
the Rashba SOC strength, we observe the scaling behavior 
Tg{E) « /3{E)Tq « j3{E)^, meaning the spin relaxes af¬ 
ter a hnite number of precession periods P (P ^ 4.5 close to 


the Dirac point), see Supplementary material. This suggests 
that dephasing is the limiting factor of spin lifetimes in the 
ultraclean case. We hnally note that by taking = 5 /reV 
(electric field of 1 V/nm jU), a spin lifetime of Tg ~ 1.4 ns 
is deduced at the Dirac point, whereas at higher energies Tg 
could reach about 7 ns. 


Discussion 


Our results show a clear transition between two different 
regimes of spin relaxation, mediated solely by the scattering 
strength of the electron-hole puddles. For hBN substrates, 
spin relaxation is dominated by dephasing arising from an 
effective energy broadening induced by the puddles, and Tg 
scales with Tp. In contrast, for Si02 substrates dephasing 
is limited by motional narrowing, leading to a DP regime 
with Tg cx 1/Tp. Remarkably, both regimes exhibit similar 
values of Tg at the Dirac point and a similar M-shape energy 
dependence (Fig. 3), making it a signature of spin relaxation 
in graphene for all puddle strengths. The crossover between 
both mechanisms occurs when Tp ~ Tq, which might have 
been realized in some experiments. This could explain some 
conflicting interpretations of experimental data in terms of 
either Elliot-Yafet or Dyakonov-perel mechanisms in. 

Our findings suggest alternative options for determining 
the spin relaxation mechanism in graphene from experimen¬ 
tal measurements. Indeed, the typical approach, to examine 
how Tp and Tg scale with electron density and to assign either 
the EY or DP mechanism accordingly, is not always appro¬ 
priate. Eor example, the EY mechanism in graphene is given 
by Tg (X E^ ■ Tp, such that Tg and Tp would scale oppositely 
with respect to electron density if Tp cx 1/Ep 0. Similarly, 
for our results the scaling of Tp and Tg with energy suggest 
an EY mechanism near the Dirac point and a DP mechanism 
at higher energies, but Pigs. 3 and 4 indicate a richer behav¬ 
ior. Therefore, to determine the spin relaxation mechanism 
it would be more appropriate to study how Tg and Tp scale 
with defect density or mobility at each value of the electron 
density. 

Pinally it should be noted that our simulations are per¬ 
formed using a constant Rashba spin-orbit coupling interac¬ 
tion, Afl, which is different from the experimental situation 
where Xp will be increased at higher charge density owing 
to larger applied external electric held. This might explain 
why, especially for hBN substrate, the simulations show a 
larger variation of Tg in energy than the gate voltage depen¬ 
dent spin lifetimes reported in the experiments Ham. 
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METHODS 

Model of homogeneous SOC and electron-hole puddles. 

The tight-binding (TB) Hamiltonian for describing spin 
dynamics in graphene is given by 

■H = -'yo'^cjcj +i^Vi c|s • {dkj x d^k)cj 
(iJ) 

+ • (s X (2) 

{i-i) 

where 70 is the nearest-neighbor 7 r-orbital hopping, Vj is the 
intrinsic SOC, and Vr is the Rashba SOC. In the low-energy 
limit, this Hamiltonian is often approximated by a contin¬ 
uum model describing massless Dirac fermions, = 

hvp^ ■ k + Xjd^Sz + Xr{<t X s)^, where vp is the Fermi 
velocity, Uk is the momentum, s{a) are the spin (pseu¬ 
dospin) Pauli matrices, Xr = |Vr, and A/ = 3-\/3V/. 
The value A/ = 12 peV is commonly used for the intrin¬ 
sic SOC of graphene i) while the Rashba SOC is electric 
field-dependent. Here, we let Xr = 37.4 p,eV, taken from an 
extended sp-band TB model for graphene under an electric 
field of a few V/nm ^ |5l. Higher-order SOC terms in the 
continuum model beyond allow an extension to higher 
energy f?!). 

Spin dynamics methodology. 

The time-dependent spin polarization of propagating 
wavepackets is computed through ll35]l 

p(F I) {'^it)WSiE-n) + s{E-n)s\'^{t)) 

^ ’ ’ 2{'i>{t)\S{E-n)\'i>{t)) ’ 

where s are the Pauli spin matrices and S{E — H) is the 
spectral measure operator. The wavepacket dynamics are 
obtained by solving the time-dependent Schrodinger equa¬ 
tion ioi, starting from a state = 0)) which may have 
either out-of-plane (^-direction) or in-plane spin polariza¬ 
tion. An energy broadening p is introduced for expand¬ 
ing 5{E — H) through a continued fraction expansion of 
the Green’s function |0Ol, and mimics an effective disorder. 
This method has been used to investigate spin relaxation in 
gold-decorated graphene lf35l . Here, we focus on the expec¬ 
tation value of the spin z-component P^{E,t) = P±{E,t) 
and the spin x-component Px{E, t) = Py {E, t). 
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